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Abstract 

We consider cyclic Lotka-Volterra models with three and four strategies where at every inter¬ 
action agents play a strategy using a time-dependent probability distribution. Agents learn from 
a loss by reducing the probability to play a losing strategy at the next interaction. For that, an 
agent is described as an urn containing (5 balls of three respectively four types where after a loss 
one of the balls corresponding to the losing strategy is replaced by a ball representing the winning 
strategy. Using both mean-field rate equations and numerical simulations, we investigate a range of 
quantities that allow us to characterize the properties of these cyclic models with time-dependent 
probability distributions. For the three-strategy case in a spatial setting we observe a transition 
from neutrally stable to stable when changing the level of discretization of the probability distri¬ 
bution. For large values of ( 3 , yielding a good approximation to a continuous distribution, spatially 
synchronized temporal oscillations dominate the system. For the four-strategy game the system is 
always neutrally stable, but different regimes emerge, depending on the size of the system and the 
level of discretization. 

PACS numbers: 05.70.Ln,05.50.+q,64.60.De,89.75.Kd 
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I. INTRODUCTION 


Evolutionary game theory is nowadays a well established approach to model and study 
complex biological and ecological systems dm. Of special importance are thereby spatial 
systems, as they give rise to novel and rich phenomena, ranging from the formation of com¬ 
plicated space-time patterns, in the form of spirals or cluster coarsening, to the emergence 
of nested ecological niches [5H8]. Most studies published in this context focused on spatial 
cyclic games with three [9H32] or four species [331(38] , but there have also been some attempts 
to generalize this to more complex networks with a larger number of species [7J ESI E2HSZ] • 
In general, these competition models assume that an agent will play a single strategy when 
confronted by another agent. The act of playing a single strategy is known in classical game 
theory as a pure strategy [58] . In the literature on population and evolutionary dynamics a 
strategy is associated with a single species and/or genotype. Depending upon the model dy¬ 
namics a losing agent will either be removed completely or it will replace the losing strategy 
with the winning pure strategy. 

In game theory, other than pure strategies, there is also the notion of mixed strategies 
[58] ■ A mixed strategy is when at every interaction an agent picks and plays one of the 
possible pure strategies using a probability distribution. Mixed strategies are sometimes 
encountered in nature [59, '60j, but can be seen most readily in social systems where decision 
making is important [HI.) ;®] • Economics is another area where game theory with either pure 
or mixed strategies has been studied. A series of economics papers focus on modifications 
of the three-species rock-paper-scissors model (see, for example, [631(67] ). whereas others tie 
three-strategy cyclic domination to the Public Goods game [681173] , Some papers discuss 
spatial games in an economics setting, but usually only games with two pure strategies 
(Prisoner’s Dilemma, Snowdrift, Hawk and Dove) are considered [7JH77]. 

In this paper we study versions of the cyclic three- and four-strategy games where agents 
are selecting a strategy from a time-dependent probability distribution. Learning from a 
recent loss, an agent changes their probability distribution in such a way that it becomes 
less likely to play the losing strategy at the next interaction. As the emerging space-time 
properties in a spatial setting are of special interest, we focus in this work on agents that 
live on a one-dimensional lattice and only interact with their two nearest neighbors. Some 
results are also presented for the well-mixed situation without an underlying spatial system. 
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In our implementation we describe every agent as an urn that contains f3 balls of three 
respectively four different types, corresponding to the different pure strategies. At every 
interaction one of these f3 balls is chosen randomly and the related strategy is played. If 
that strategy loses, the ball is replaced by a different ball representing the winning strategy. 
In that way the losing agent is less likely to play at the next interaction the losing strategy 
again. It should be noted that the value of 3 is a measure of discretization of the probability 
distribution, with the limit j3 —)• oo yielding a continuous distribution. As we discuss in 
the following, changing the value of /3 has a strong impact on the properties of systems with 
cyclic domination. For example, for the three-species strategies we observe a transition from 
a neutrally stable system to a stable system when increasing /3. In addition, in the limit of 
1 the system synchronizes, yielding spatially extended coherent temporal waves. When 
considering four strategies, one always has a neutrally stable system, but the average time 
for strategy extinction displays different regimes, depending on the value of f3 and on the 
system size. 

Our paper is organized in the following way. In the next section we present our models 
in more detail. Sections III and IV are devoted to our results, first for the three-strategy 
case and then for the four-strategy model. In order to elucidate the properties of our 
mixed-strategy spatial systems with cyclic domination and a time-dependent probability 
distribution we study a range of quantities: time-dependent densities, space-time covariances 
and related length scales, and average times for strategy extinction. We summarize and 
conclude in Section V. 

II. MODELS 

In the symmetric Lotka-Volterra model reactions are taking place with rate A in a cyclic 
way between M different species j?8j|79]: 

W + x i+1 4 Xi + Xi ( 1 ) 

where X t is an agent of species i, with i = 1, • • • , M and the cyclic identification M + 1 = 1. 
In the language of population dynamics, every species has one prey and is the prey of a 
single other species. The three-species case is special, as a reaction takes place whenever 
two agents of different types are selected. For cases with more than three species one has 
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mutually neutral species that do not interact [79| 80]. Cyclic Lotka-Volterra games can be 
played in the well-mixed case as well as on lattices where a variety of situations have been 
considered. For example, one can impose a strict site restriction with exactly one agent 
occupying each lattice site and reactions taking place between neighboring sites pi 
This can be coupled with exchanges of agents sitting on nearest neighbor sites, in order to 
provide a mechanism for mobility mm- Sometimes the site restriction is dropped so that 
every site can hold a variable number of agents that can diffuse and interact with agents on 
neighboring sites [IS]. In the related systems with May-Leonard dynamics, the reaction (|TJ) 
is replaced by two separate reactions P, ECE] (usually with different rates) 


Xi + X i+1 —> Xi + 0 


Xi T 0 y Xi T Xi 


( 2 ) 


where 0 indicates an empty site. In this scheme the number of agents in the system is a 
stochastic quantity that is only conserved on average. 

We consider in the following versions of the three- and four-strategy cyclic Lotka-Volterra 
game with mixed strategies where at every interaction agents play one of the possible strate¬ 
gies using a time-dependent probability distribution. Imagine an agent as an urn that con¬ 
tains /3 balls of three (for the three-strategy version) or four (for the four-strategy case) 
different types. For a well-mixed system without a lattice structure we first choose two 
agents/urns at random out of a total of N agents/urns before choosing randomly a ball out 
of each selected urn. Depending on the types of balls selected, a reaction may take place 
with rate A = 1 following the scheme (JTj) . If the strategy played by one of the agents is 
beaten, then the losing ball is replaced by a ball of the winning type before the balls are 
put back into the urns. The losing agent therefore changes the probability distribution as 
a result of the loss by increasing (decreasing) the probability to play the winning (losing) 
strategy at the next interaction. For a spatial game we only consider the case with exactly 
one agent at every lattice site. In order to start an interaction we select an agent and one of 
their neighbors at random and then proceed as for the well-mixed case. In our simulations 
we define one time step to be N/3 proposed updates. 

One can view this scheme in a spatial setting as a version of the spatial cyclic Lotka- 
Volterra model with multiple occupancy of each site, but there are important differences. In 
the model discussed in [16] in the context of population dynamics, individuals not only inter- 


4 


act with agents on other lattice sites, they also diffuse by jumping with a given probability 
to one of the neighboring sites. As a result the number of individuals at each site fluctuates 
and only the average spatial density is constant. In our case, only interactions take place 
so that the number of balls at each site (which corresponds to the number of individuals at 
a given site in the model of [16]) is strictly conserved. As we will see in the following, this 
additional conservation law has a huge impact on the properties of our system. 


III. THE MIXING OF THREE STRATEGIES 

As already mentioned, the situation with three strategies is rather special, as every time 
two different strategies are played, there will be a losing and a winning strategy. In the 
following we show, both in the well-mixed case as well as on the ring, that for continuous 
probability distributions a synchronization of the strategies played by the different agents 
takes place. For small number of balls per agent a transition in the stability properties of 
the lattice system is observed. 

A. Mean-field equations for the well-mixed case 

We first consider the mixed three-strategies game in a well-mixed system without spatial 
structure. Agent j is characterized by the number of balls of each type at their possession: 
(aj, bj , Cj ), with dj + bj + Cj = /3, where aj is the number of balls of type A. The state of the 
system of N agents is then given by the number of balls of each type in possession of each 
agent, i.e. by the configuration 

{ b n7 C n ) j* {(otl, bi, Cl), - • ■ , ( (lj , bj , Cj ), ■ ■ ■ , ((X 7 V 5 bjy , Cjy) } . (3) 

The interaction scheme (|T|) rates A = 1 directly translates into the following Master 
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equation for the probability that the system is in state {(a n , b n , c n )} at time r + 1. 


(fl n , 6 n , C n ) } , T -)- 1 ) 

1 2 
]PN(N- 1) 


+ a *( 6 i + 1 ) jP (''' > ( a *> Ci)> ■ • • > ( a i - !, bj + 1 , Cj), • • • ; t) 

+ X] S 6 <( c i + 1 ) p (''' > ( a *> 6 i» c *)> • • • > ( a i> b j “ !, 9 + !), • ■ • ; T ) 

+ 5ZS Ci ( a J + 1 ) p ('" >( a i>bi,Ci),--‘ ,(aj + l,bj,Cj-l),--- ;t) 

I 2 _ A _ 

+[! - - 1) ^ + 6<c i + c ‘ a J')]- P ({( a "> c ™)} 5 r ) 

( 4 ) 

where the factor ^ takes into account the number of combined choices of the two agents 
i and j. The hrst three sums represent the reactions 0 through which the system enters 
the state {(a n ,b n ,c n )}, while the last term represents the case where no reaction happens 
(when balls of the same type are pulled out for both agents) and the system stays in the 
same state. 

We define the density of ball type A for agent j at time r to be 


Aj(r) = °^ p {{{ a n,b n ,c n )} ;r) (5) 

{(&n ibn ?Cn)} 

with similar equations for ball types B and C. Neglecting correlations, Eqs. ® and Q 
yield the equation (with similar equations found for the other ball types through symmetry) 


A k {r + 1) - A k (r) = ^ ^ (P({(a n ,b n ,c n )}-,r + 1) - P({(a n ,b n ,c n )} ]r)) 

{ (&n ,&n ? c n)} 

_ l 2 _ A 

= 5Z 03 N ( N - 1) J2( aibk ^ C i a k)P{{{ a m b n , C n )} ; t) 

{(&7i , c n)} i^k 


1 

iV/3 



B k (r) - 2 



A k {r) 

( 6 ) 


Letting t — jh and taking the continuum limit /3 —> oo so that A k {r + 1) — A k (r) —$■ 
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-j^-pd t Ak(t) finally leads to the mean-field equations 


( 7 ) 


d t A k (t) = 2 ( ) B k {t)~ 2 ( ) Ak $) 

d,B t (t) = 2 ( ) C k (t) - 2 ( J2 A ‘(t) ) B„(t) 

a,C k (t) = 2 ( —]T Q(t)) 4t(t) - 2 ( -i- ]T Bi(t)) c„(t) 

V i¥=k / \ i^k / 

As we show in Appendix A, the usual mean-field rate equations for the well-mixed rock- 
paper-scissors model are recovered if one takes in addition the limit of infinitely many agents, 
N —> oo. 

Let us have a closer look at the simple case of two agents N = 2. Eqs. (JT|) can be easily 
written as 

dt,A\ = 2(A 2 i?i — C 2 A 1 ) , d t A 2 = 2 (AiB 2 — CiA 2 ) 

M = 2(B 2 C x - A^) , d t B 2 = 2[B X C 2 - AM (8) 

dtC\ = 2(C 2 A\ - BM , d t C 2 = 2(C' 1 A 2 - BM 

From d t (A] — A 2 ) we hnd 


dt{A\ — A 2 ) — 2 [( A 2 Bi — C 2 Ai) — ( AiB 2 — CiA 2 )\ — —2(A 1 — A 2 ) 
as Ai + Bi + Ci = 1. Solving this differential equation yields 

Ai(t) - A 2 (t ) = (ai - a 2 )e~ 2t 


(9) 


( 10 ) 


where a\ = Ai(0) and a 2 = A 2 (0) are initial conditions. Similar equations are obtained for 
ball types B and C. It follows that whatever the initial conditions the difference in ball 
densities of the two agents vanish exponentially fast, yielding a synchronization of the two 
agents. Although not presented here, the synchronization of more than two agents has also 
been studied in some detail, see [82]. 

We can use the previous results to eliminate the densities of agent 2 and end up with the 
following three coupled non-autonomous differential equations for the ball densities of agent 
1 : 

9tA\ = 2(A 1 B 1 — C\Ai) + 2 [(a 2 — a\)Bi — (c 2 — c\)A{\ e 2t 

d t B 1 = 2{B l C 1 - AM + 2 [( b 2 - ln)C x - (a 2 - a^B,} e~ 2t (H) 

d t C\ = 2(C' 1 A 1 - BM + 2 [(02 - Cl )Ai - (b 2 - bJCi] ^ 2t 
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We note that in the long time limit t oo we recover the mean held rate equations for the 
well-mixed rock-paper-scissors game. 



FIG. 1: (Color online) Time-dependent densities for a system of two agents obtained through 


numerical integration of Eqs. (11) (full black line: agent 1, green dashed line: agent 2) and 
compared with results from stochastic simulations with (3 = 10' (red circles: agent 1, blue squares: 
agent 2). The same initial conditions Ai(0) = 0.8, -Bi(O) = 0.1, ^(0) = 0.1, and B 2 (Q) = 0.8 
are used for both methods. In (d) the symbols represent the absolute values of the difference 
between strategy densities in the simulation data, while the lines are obtained from Eqs. An 

exponentially fast synchronization is observed. The dot-dashed blue line indicates an exponential 
of slope —2. 


In Fig. [I] we compare these theoretical results (lines in the figure), obtained through 
numerical integration of the Eqs. GD- with results from stochastic simulations (symbols in 
the figure) for initial conditions kLi(0) = 0.8, £q(0) = 0.1, A 2 { 0) = 0.1, and £> 2 ( 0 ) = 0.8. 
We find that the numerical integration results agree with the stochastic simulation very 
well; in analogy with the three-species well-mixed model J53J we expect the fluctuations 
for an individual agent to be proportional to /3~ 1 ^ 2 . Even though the initial conditions are 
very different, the synchronization between the two agents is very rapid and the differences 











between the corresponding densities vanish exponentially, see Fig. p. 

B. Numerical simulations on the ring 

The simplest spatial system is the one-dimensional lattice with periodic boundary con¬ 
ditions. As discussed previously, every site of the ring is occupied by one agent who has at 
their disposal f3 balls. In the initial state every ball is assigned with the same probability 
one of three possible types corresponding to the three strategies A, B , and C. Once the 
system has been prepared in that way, pairs of neighboring sites are randomly selected and 
interactions take place following the scheme described above. 

The dynamics can be readily visualized through space-time plots as those shown in Fig. 
[2j Inspection of these plots for various numbers of balls (3 reveals an interesting transition 
in the shape of the space-time patterns. Whereas for three balls or less the system behaves 
like the standard Lotka-Volterra rock-paper-scissors game on the ring with immobile agents 
(which corresponds to f3 = 1) and exhibits coarsening processes that end when only one 
strategy fills the complete lattice (see Fig. |2^i), for (3 > 4 a tiling structure appears where a 
tile indicates that a part of the system is dominated by one of the three strategies for a finite 
amount of time, see Fig. !>• It follows that every agent changes the most likely strategy 
after some time and that it becomes difficult for one strategy to dominate the system. This 
tiling structure is reminiscent of very similar patterns that are encountered when allowing 
in the one-dimensional rock-paper-scissors model for swapping of particles as an efficient 
mechanism for mobility [13]. As discussed further below, see Fig. [3j the tiling structure 
promotes coexistence and stabilizes the system, with only rare large stochastic fluctuations 
causing the finite system to go to an absorbing fixed point. When we further increase (3, 
the tiles, corresponding to regions where one strategy dominates locally, decrease in size, see 
Fig. [2]: for the case with (3 = 6. This is accompanied by the emergence of grayish patches 
that indicate spatial regions where in the probability distribution the three strategies have 
similar weights. Finally for large values of /3, corresponding to a probability distribution 
that approximates a continuous distribution, the system rapidly synchronizes and spatial 
extended coherent temporal waves are formed, as shown in Fig. [2jl for the example of 
/3 = 100 balls. 

In order to relate the coarsening and tiling space-time patterns to different extinction 
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FIG. 2: (Color online) Space-time plots for three-strategies games on a ring with different numbers 
of balls (3 on each site: (a) (3 = 3, (b) /3 = 4, (c) (3 = 6, and (d) f3 = 100. Time progresses from top 
to bottom. For (a)-(c) 1000 time steps are shown for a system composed of 1000 sites, whereas for 
(d) we show 500 time steps for a system with 500 lattice points. To determine the color of a lattice 
site we use the RGB color model and map the percentage of A, B, and C to the percentage of the 
colors Red, Green and Blue respectively. In this scheme a site that contains the same number of 
balls for all three species is assigned a grayish color. 
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FIG. 3: (Color online) Average lattice domination time T as a function of the total number of 
balls N/3 in the system. When changing the number of balls (3 per agent, a transition takes place 
between an algebraic dependence on the total number of balls in the system for f3 < 3, see panel 
(a), indicating a neutrally stable system, and an exponential dependence for (3 > 4, characterizing a 
stable system, see panel (b). Each data point results from an average over 2000 runs with different 
realizations of the noise. Error bars are smaller than the symbol sizes. 

regimes, we measure the average lattice domination time T, i.e. the average time at which 
only one of the strategies remains in the system. Fig. [3] reveals that the average lattice 
domination time changes its dependence on the system size at the transition gleaned from 
the space-time plots. For (3 < 3 the average lattice domination time increases algebraically 
with the total number of balls, T ~ (N (3) b , as shown by the straight lines in the log-log plot 
in panel (a) of Fig. [3j with the exponent b = 1.93(2), 1.90(2), and 1.77(2) for j3 = 1, 2, 
and 3, respectively. As shown in Ref. [9j EH E5|, there is a direct correspondence between 
the system size dependence (exponential, algebraic, or logarithmic) of the lattice domination 
time of a system and its stability properties (stable coexistence, neutrally stable coexistence, 
or unstable coexistence). The observed algebraic dependence therefore indicates that the 
system is neutrally stable. We also note that the lattice domination time for a fixed value of 
N/3 decreases for increasing (3, indicating that the system becomes less stable. This trend 
is reversed when we enter the tiling regime, see Fig. if. as now for fixed N/3 the lattice 
domination time strongly increases with (3. In fact, the dependence of T on N(3 changes to 
an exponential dependence for (3 > 4, as the system becomes a stable system [83] . 

Another quantity that allows us to gain insights into this transition is the total space-time 
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FIG. 4: (Color online) Spatial covariance at time t = 1000 for different numbers of balls /3 
(increasing from 1 to 6 from (a) to (f)) and different system sizes. The data result from averaging 
over half a million independent runs. 

covariance 


n , C A a(x, t ) + C B b(x, t) + C C c(x, t) 

L'o(2b t) — g 

(12) 

C A A(x,t ) = y ^Ai(t)A i+x (t) - n A (t)ii A (t) 

(13) 


i 


and similar expressions for CsB(x,t) and Ccc(x,t). Here HA(t) = y? Yi A(£) and similarly 
for n B and /z<> Fig- [4] shows this quantity after 1000 time steps since the preparation of the 
system for f3 ranging from 1 to 6 and system sizes between N = 500 and N = 10000. For 
/3 < 3 the covariance displays the expected behavior for systems with domain ordering, with 
strong finite size effects for small systems, as in many instances runs have already reached 
their final state (with one strategy filling the whole system) at t — 1000, and an exponential 
decay with the distance x for larger sizes when the coarsening system is still far from its 
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final state. After the transition to the tiling regime, see (d) to (f) in Fig. |4j the initial decay 
becomes system size independent and is followed by a shoulder (see the data for [3 — 6 in 
Kg- §)• We relate this behavior to the typical sizes of the tiles whose spatial extensions 
are rather small and decrease with an increase of [3. 

C. Spatial mean-field equations and synchronization 

As we saw in Fig. [2]i synchronization in space coupled with temporal oscillations sets in 
for large number of balls (3 . As in the limit (3 —» oo the probability distribution becomes 
continuous, we can capture this effect through spatial mean-held equations. 

Our starting point are the mean-held equations ([7]) for N agents in the well-mixed case 
with [3 —> oo that need to be adapted to the spatial setting of a one-dimensional lattice 
where the agent on lattice site k interacts exclusively with the agents located on the neigh¬ 
boring sites k — 1 and k + 1. Each sum in ([7]) over i ^ k then reduces to two terms, so that 
we obtain for A k (t) (with similar equations for B k (t) and Ck(t )): 

dtA k = (A fc+ i + Ak-i)Bk — (Ck+i + Ck~i)Ak . (14) 

Introducing the finite difference A 2 A k = ( A k+ \ — A k ) — ( A k — A k _ 1 ) allows to recast this 
equation in the form 

dtA k = (A 2 A k )B k — (A 2 C k )A k + 2 (B k A k — C k A k ) . (15) 

Taking the spatial continuum limit, where we approximate A by ad x , yields the following 
set of partial differential equations 

d t A(x , t ) = ( d 2 A(x , t )) B(x, t ) — (<9 2 C(;r, t )) A(x, t) + 2 (B(x, t)A(x , t) — C(x , t)A(x, t )) 

d t B(x, t ) = ( dlB(x , t )) C(x, t ) — [dlA(x : t )) B(x , t) + 2 ( C{x , t)B(x, t ) — A(x, t)B(x , t)) 

d t C(x,t ) = (d 2 x C(x,t)) A(x,t ) — ( d 2 B(x,t )) C(x,t) + 2 (A(x,t)C(x,t) — B(x,t)C(x,t)) 

(16) 

where we set the length scale a — 1. This set of equations can straightforwardly be general¬ 
ized to d dimensions. The terms d 2 describe the diffusion of strategies through interactions, 
as the loser of an interaction changes their probability distribution in favor of the strategy 
against which they lost. The remaining terms are nothing else than the mean-held equations 
of the normal three-species cyclic game in the well-mixed case. 
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FIG. 5: (Color online) Space-time plots from (a) numerical integration of the spatial mean-field 


equations (16) and (b) the numerical simulation of a system of 100 lattice sites over 100 time steps 
with f3 = 10 6 balls at each site. In both cases the same initial condition was used. Time 
increases in the downward direction. The system rapidly synchronizes and coherent waves are 
formed. 

In Fig. [5^i we show the numerical solutions of this set of equations for the initial condition 


A{x,t = 0) = -Q[N /5 - |a: - N/2\\ 
B(x,t = 0) = le[Af/4 - |ap - JV/2|] 


(U) 


C(x,t = 0) = 1 — A(x, t — 0)■*- B(x, t — 0) 

where ©[---] is the Heaviside step function. In this initial state large segments of the 
system are occupied by agents that play the same initial strategy. For comparison we show 
in Fig. [ 5)0 a numerical simulation for (3 = 10 6 balls and the same initial condition. As 
expected a system with such a large number of balls is well described by the mean-held 
equations. The space-time plots in Fig. [5] show that the spatial system with a continuous 
probability distribution synchronizes very rapidly, yielding spatially extended regions where 
the probability distribution of every agent is very similar. As a result the strategies in the 
whole system coherently oscillate in time. This synchronization effect, which is independent 


of the initial condition, is readily understood from Eqs. (16) by remarking that the diffusion 


terms efficiently smooth out the spatial inhomogeneities in the probability distributions until 
only the well-mixed terms describing the standard rock-paper-scissors interactions matter. 
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IV. THE MIXING OF FOUR STRATEGIES 


It results from the reactions (|TJ) that in the cyclic Lotka-Volterra scheme with four strate¬ 
gies, pairs of mutually neutral strategies are encountered, as the strategies A and C (B 
and D) do not compete against each other [803 • I 11 the case of pure strategies, this part¬ 
nership formation yields in a spatial game the formation of domains composed of neutral 
partners [33] as an agent with a given strategy takes advantage of the fact that its neutral 
partner plays a strategy that beats the strategy against which the agent would lose. This 
guiding principle also holds true when considering a four-strategies mixed game with time- 
dependent probability distributions. Still, remarkable changes in the domain structure and 
in the mean lattice domination time take place when changing the level of discretization of 
the probability distribution by increasing /3. 

A first impression of the changes that happen when the number of balls /3 is increased 
can be gained from the space-time plots shown in Fig. [6] For /3 = 1 we see the formation 
of regions dominated by neutral pairs (red and blue vs green and black). As there is no 
mechanism for mobility, red and blue (green and black) single-species regions become stuck 
within one another, forming superdomains of neutral partners [7E1 E9]. As a result a winning 
strategy can invade the region of a losing strategy only until it hits a patch occupied by the 
partner strategy of the losing strategy. This results in the zig-zag like structures where the 
different strategies (in the order red, black, blue, green) dominate one after the other over 
a region of the lattice. These regions grow in extent after each change of strategy, yielding 
ultimately a lattice occupied by only one of the partnerships (either red and blue or green 
and black). As we increase /3, domains of neutral pairs become effectively mixed, and a third 
strategy has a higher probability to invade a superdomain occupied by a given alliance. This 
results in a much slower growth of the neutral pair domains. Further increasing f3 causes 
the neutral species pairs to become very well mixed, resulting in two types of neutral-species 
domains that overall look purple and dark green. These neutral-species domains continue 
to compete against each other, yielding a slow coarsening process. Within these coarsening 
domains cyclic processes continue unabated, as revealed by the appearance of localized wave 
patterns for very large f3, see Fig. [6ji. These waves are very much like those encountered in 
the three-species case, albeit of smaller extent, only that now one cycles through four types 
of balls instead of three. 
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FIG. 6: (Color online) Space-time plots for four-strategies games on a ring with different numbers 
of balls j3 on each site: (a) /3 = 1, (b) /? = 10, (c) /3 = 20, and (d) /3 = 50. Time progresses from 
top to bottom. 1000 time steps are shown for systems composed of 1000 sites. To determine the 
color of a lattice site we use the RGB color model and map the percentage of A, B, and C to the 
percentage of the colors Red, Green and Blue respectively. This is adequate as the percentage of 
D is readily obtained from the fact that the sum over all four densities is 1. 


For a more quantitative discussion we turn again to the space-time covariance of the form 


1 


Cxr(x,t) = — ^ X,(t)Y, +x (t.) ~/t x (t)fiy(t) 


(18) 


which yields the self-covariance if the species X and Y are the same (see Eq. (13)) and 
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the cross-covariance otherwise. We then define an individual spatial covariance through the 
equation 

Ci{x,t) = ^ (C AA (x,t) +C BB (x,t) + C C c(x,t) +C DD (x,t)) . (19) 

A different spatial covariance can be obtained when we do not distinguish anymore between 
neutral strategies but consider them to form a unique group, i.e. strategies A and C together 
form the group denoted by A, whereas B and D make up the group B. We call the space-time 
covariance for these larger groups 

C n {x, t) = ^ (C AA (x, t ) + C BB (x, t)) (20) 

the neutral spatial covariance. From these two quantities we can extract two different time- 
dependent length scales, Lj(i) and L n (t), that provide some insights into the ordering of the 
superdomains formed by neutral partners. These lengths are obtained from the intersection 
of the normalized covariance, Ci(x,t)/Ci(0,t) and C n (x,t)/C n (0,t), with a line of a constant 
value k, i.e. Ci(Li(t),t)/Ci(0,t) — k and similarly for L n (t). We use in the following k = 0.5 
after carefully checking that the qualitative features discussed below do not depend on the 
precise value of k. 

Fig. 0 shows the time evolution of these two lengths for different values of (3. Let us 
first look at the pure case with /3 = 1, as we can compare for this case our lengths with 
those discussed previously in the literature. As shown by Frachebourg et al ng one needs 
to consider two different types of domains for the four-species Lotka-Volterra model with 
immobile particles on a one-dimensional chain: the pure domains occupied by a single species 
and the superdomains shared by two non-interacting species. Starting from a fully disordered 
state, the pure domains increase as t 1//3 , whereas the superdomain size, i.e. the distance 
between active interfaces, is proportional to t 2 / 3 . As seen in panel (a), for the case [3 — 1 
both lengths Lj(f) and L n (t) provide essentially the same information on the superdomains: 
they are proportional and both display an algebraic increase with an exponent 2/3 before 
saturating due to finite-size effects. 

Increasing f3 yields a slowing down of the coarsening process as manifested by a gradual 
decrease of the exponent governing the growth of the two lengths. For /3 = 5, for example, 
the effective exponent before the transition to the saturation regime is close to 0.45, see 
panel (b). This decrease of the exponent describing the long-time regime before saturation 
continues when further increasing (3, with a value of 0.14 for (3 = 15 and 0.05 for (3 = 30, 
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FIG. 7: (Color online) Time-dependent lengths extracted from the space-time covariance for the 
four-strategies mixed model on a one-dimensional lattice. The different panels show results for 
different values of /3: (a) /3 = 1, (b) /3 = 5, (c) /3 = 15, (d) (3 = 30. Li respectively L n is obtained 
from the individual respectively neutral spatial covariance. In (a) the system size is N = 5000, 
whereas in all other panels data for a system of 2500 sites are shown. The data result from averaging 
over typically a few thousand independent runs. In (a) the dashed line indicates an exponent of 
2/3, whereas in (c) and (d) the short red lines indicate an algebraic growth with an exponent of 
1 / 2 . 

see panels (c) and (d). This very slow increase for large f3 is not size-dependent and should 
therefore not be confused with the finite-size plateau seen for example in (a) for (3 = 1. 
The very weak increase of the correlation length for large /3 and large t instead indicates 
that domain growth almost comes to a standstill due to strong mixing effects. Interestingly, 
for large /3 this long-term regime is preceded by another regime where L n increases as a 
square-root of time (indicated by the red segments in (c) and (d)), whereas Lj displays 
some non-trivial features that reflect the complicated ordering processes seen in the space- 
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time plots. We tentatively identify this regime with the initial formation and growth of the 
neutral-species superdomains followed by a coarsening process where essentially only two 
types of domains compete against each other. 
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FIG. 8: (Color online) Average lattice domination time as a function of the total number of balls 
N (3 in the system. Each data point results from an average over 2000 runs, the standard error 
being smaller than the sizes of the symbols. 

Another way to characterize these systems is through the study of lattice domination 
events. The lattice domination time shown in Fig. [8] is the mean number of time steps 
needed until only one neutral species pair remains in the system. Interestingly, different 
regimes also show up in the lattice domination time when changing the value of (3. For 
(3 = 1, i.e. the case of pure strategies, the lattice domination time increases with the system 
size as N 1A0 |37j . making this a neutrally stable system. Increasing (3 slightly changes the 
exponent (for (3 = 4 its value is 1.68), but has no other effect on the lattice domination time. 
For (3 > 5, however, different regimes emerge as a function of the system size, as shown 
in Fig. [H£i. Plotting the lattice domination time against the total number of balls in the 
system, N(3, reveals a first regime where the lattice domination time increases algebraically 
with an exponent close to 2. This is followed by a crossover to a second algebraic regime 
with a much larger exponent (for example, for /5 = 30 the exponent is 5.60). For (3 not 
too large, this crossover takes place at rather similar values of N(3, i.e. the larger (3 is, the 
smaller the needed system size is to enter the second algebraic regime. Fig. [8]d displays 
another change of behavior for values of (3 > 20, as the crossover is then shifted to larger 
values of Nfi when increasing (3. For the smaller sizes, i.e. before the crossover, the data for 
different (3 values fall on one common curve with an exponent of approximately 1.5. 
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The emergence of two regimes in the lattice domination time for /3 > 5 can be related to 
different types of extinction events linked to the prevailing domain structure. On the one 
hand, in small systems extinction can take place at rather short times since the preparation 
of the system, due to the formation of a few superdomains followed by a coarsening of 
these domains. For larger N many such domains are formed, resulting in complicated 
processes dominated by triangular (zig-zag) space-time patterns as those seen in Fig. & 
for — 10. This periodic cycling through triangles of all four ball types yields a very slow 
coarsening process. A further increase of f3 results in the replacement in this second regime 
of the triangular structures by very long-lived wave patterns, due to some spatially localized 
synchronization, that dominate the purplish and dark green domains in Fig. [6]T. Domains 
grow very slowly in that regime, see the late-time behavior of the growth length shown in 
Fig. [ 7 JI, and extinction events are only encountered at very late times. 

In [37] we showed for the four-species model with pure strategies ( f3 = 1) that much 
can be learned about extinction events when studying the probability distribution of the 
lattice domination time. We expect this to be also true for the more complicated cases with 
(3 large. However, the extremely large values of the domination times make it impossible 
with the resources at our disposal to perform enough runs for a reliable computation of this 
probability distribution. 


V. CONCLUSION 

Systems composed of multiple species that interact in a cyclic way have been at the center 
of a multitude of studies in recent years. Mostly discussed in the context of evolutionary 
game theory and population dynamics, these systems allow to understand some of the generic 
properties arising from non-trivial interactions in ecological systems. 

In this paper we discussed a version of the three- and four-species Lotka-Volterra model 
where the agents are using a mixed strategy, i.e. agents play a pure strategy using a proba¬ 
bility distribution every time they interact. Taking into account that agents, both in ecology 
and economics, tend to learn from past experience, we consider time-dependent probability 
distributions where a losing agent decreases the probability to play a losing strategy at the 
next interaction. In order to do so, we treat an agent as an urn containing f3 balls of three 
respectively four types, where each type corresponds to one of the three respectively four 
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strategies. If a strategy loses, a ball of the corresponding type is replaced by a ball with 
the winning strategy. The number of balls (3 in the urn therefore measures the level of 
discretization of the probability distribution, with the limit (3 —> oo corresponding to a 
continuous distribution. 

As our study revealed, some remarkable changes take place when changing the level 
of discretization. For the three-strategy case on a one-dimensional lattice we observe a 
transition between neutrally stable for (3 < 3, where the average time needed for one strategy 
to completely pervade the system grows algebraically with the system size due to coarsening 
of single domains, to stable for /3 > 3, where this time increases exponentially with the 
extent of the system, thus indicating the existence of a stable attractor in the coexistence 
region. This transition gives rise to a change of space-time patterns and the emergence of a 
tiling structure where strategies dominate for a finite amount of time over certain regions of 
the lattice. In the limit (3 1, when the probability distribution approximates a continuous 

distribution, this tiling structure evolves into spatially extended waves where the dominating 
strategies changes periodically. Synchronization throughout the system is also encountered 
in the case without spatial dependence and can be understood in the mean-field limit of 
an infinite number of agents using continuous probability distributions by analyzing the 
corresponding rate equations. Comparing our three strategy results to those found in [16] 
for a related multi occupancy model, we find much richer dynamics and a transition in the 
stability properties of fully occupied lattices when changing the occupancy parameter (3. 

The four-strategy case is characterized by the existence of pairs of non-interacting strate¬ 
gies. As a result agents in a spatial setting with (3 = 1 want to ally themselves with agents 
that play the complementary strategy in order to fight off the competing pair of strategies. 
A direct consequence of this rivalry between competing alliances is the formation and coars¬ 
ening of superdomains occupied by a single alliance. In contrast to the three-strategy case, 
an increase of (3 does not yield a stable system. Instead, the mean lattice domination time, 
i.e. the average time needed for one alliance to completely fill the system, always increases 
algebraically with N/3. Still, different regimes can be identified as a function of N and (3. 
For example for large values of these parameters a very slow coarsening process is observed, 
with local synchronized waves within the competing domains. 

Besides some results for the well-mixed case without spatial setting, we focus in this 
paper on the one-dimensional lattice. It is well known from cases with (3 = 1 that the 
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dimensionality of the lattice can have a huge impact on the properties of the system. Taking 
into account the already complex behavior observed in our study of the ring, we expect the 
appearance of additional intriguing features, especially for the four-strategy case, when 
considering systems with time-dependent probability distributions in two space dimensions. 

In our paper we only considered immobile agents. This is of course not a realistic descrip¬ 
tion for ecological systems. For the Lotka-Volterra type models considered here, mobility 
can be implemented through the swapping of agents occupying neighboring lattice sites. 
It would be interesting to see how the different regimes are modified when allowing for 
exchanges of agents. We expect to come back to that question in the future. 

We have restricted us to the simple three- and four-species Lotka-Volterra models with 
time-dependent probability distribution, as for these situations the properties for the case 
(3 = 1 are well understood and provide a case against which we can study changes that 
emerge when using a time-dependent probability distribution to play a strategy. Recently 
there has been a strong focus on more complicated situations, with more species and/or 
more complex interaction schemes. It is an interesting question how the properties of these 
systems change when considering a mixed strategy game with a time-dependent probability 
distribution. We expect this to be a very fruitful research avenue for the future. 
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Appendix A: Mean-field equations for the three-species case and N 


oo 


For the well-mixed case the first of the three equations in (Jt]) can be rewritten as 

a,Mt) = 2 (-^T £ Mt)) B k (t) - 2 (£ aw) m*) 

V ) \ ijtk ) 


N 


N 


N - 1 
2 

iV — 1 


Bk Ai — Ak Ci — BkAk + AkCk 

, i i 

( B k N(A) N — A k N(C)N — B k A k + A k C k ) 


(Al) 


N N A k C k — B k A k 

= 2 1 Bk W^i^ N - Ak W^i^ N + n 


where denotes a mean over all agents/urns in the system. Taking N —)• oo and 

making the index k to be continuous yields 


d t A(x, t) = 2B(x,t)(A(x,t)) x - 2A(x,t)(C(x,t)) x (A2) 


where (■ ■ • ) x = - j J indicates an average over the continuous index x. Finally, applying 


(• • • ) x on both sides of Eq. (A2) yields 


d t (A(x, t)) x = 2(B(x, t)) x (A(x, t)) x - 2(A(x, t)) x (C(x, t)) x 

d t (B(x, t)) x = 2(C{x, t)) x (B(x, t)) x - 2 (B(x, t)) x (A(x, t)) x (A3) 

d t {C(x, t)) x = 2(A{x, t)) x {C(x, t)) x - 2(C(x, t)) x (B(x, t)) x 

where the second and third equations follow from symmetry. These are exactly the mean 
field equations for the three-species well-mixed case. 
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